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Abstract 

We present a modification of Newton's metfiod to restore quadratic con- 
vergence for isolated singular solutions of polynomial systems. Our method 
is symbolic-numeric: we produce a new polynomial system which has the 
original multiple solution as a regular root. Using standard bases, a tool 
for the symbolic computation of multiplicities, we show that the number 
of deflation stages is bounded by the multiplicity of the isolated root. Our 
implementation performs well on a large class of applications. 
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1 Introduction 

Let -F(x) = be a polynomial system of N equations in n unknowns x e C". 
We are interested in x*, an isolated solution of F(x.) = 0: 

for small enough e > : { y £ C" : ||y-x*|| < e } nF"^(0) = {x*}. (1) 

Denote by A{x.) the Jacobian matrix of the system F{x.) = 0. We call x* a 
sm^uZar solution of F{x) = <^ rank(A(x*)) < n. Let m be the multiplicity of 
the isolated solution x* of F{x.) = 0. 
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Newton's method (also ealled the method of Gauss-Ncwton when N > n) 
generates a sequence of approximations for x*. If x* is nonsingular, then 
the sequence converges quadratically (i.e.: ||xfc — Xfc+i|| = 0(||xfe_i — Xfe|p)) 
to X*, which justifies its widespread usage. But otherwise, if x* is singular, the 
convergence slows down and gets lost when x^ « x*. 

The brutal force approach is to use a working precision of m x D decimal 
places to achieve D correct decimal places in the final approximation. Even as 
multiprecision arithmetic is widely available and nowadays less expensive to use, 
this approach can only work if all coefficients in the system F have their first 
mx D decimal places correct. Our goal is to restore the quadratic convergence 
of a sequence converging to an isolated singular root without imposing extra 
requirements of precision on F. 

Newton's method for singular solutions has been extensively researched. The 
research up to the mid eighties is surveyed in [18]. We classify research related 
to our work in two domains: 

1. Detection and treatment of bifurcation points. When following a so- 

lution path of a system defined by a parameter, the solution path may turn 
back or bifurcate for increasing values of the parameter. Techniques to de- 
tect and compute such bifurcation points are generally done via Liapunov- 
Schmidt reduction [1] [15]; see also [29, 30], [10], [23], and the references 
cited therein. 

2. Deflation method for polynomial systems. A symbolic deflation method 

was presented in [36], and further developed in [33], [34], and [35]. We 
discovered this approach from the reference list of [25]. 

A theoretical framework to study the complexity and numerical stability of 
Newton's method was developed by Shub and Smale, see [4], and was generalized 
to overdetermined systems in [9]. See [13, 14] for recent generalizations of this 
a-theory to multiple roots. 

The contribution of this paper is twofold. First - as announced in [40] - 
we provide a numerically stable implementation of a modified symbolic defla- 
tion method. Second, using standard bases [16], we show that the number of 
deflations needed to restore the quadratic convergence of Newton's method is 
bounded by the multiplicity. In the next section we describe our method, fol- 
lowed by an introduction to standard bases and our proof in the third section. 
Our symbolic-numeric implementations and numerical results are described in 
sections four and five. 

Acknowledgements. Our modified method with numerical results was pre- 
sented at a poster session at ISSAC'04, 6 July 2004, Santander, Spain. We are 
grateful for the reactions of Grcgoirc Lccerf - in particular for pointing at the 
work of Kunkcl - and to Erich Kaltofen and Zhonggang Zeng. Also discussions 
with Andrew Sommese and Charles Wampler at Oberwolfach helped us. 
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2 A Modified Defiation Method 



A singular root x* of a square (i.e.: N = n) system -F'(x) = satisfies 



The augmented system (2) forms the basic idea for deflation. If x* is isolated 
and corank(A(x*)) = 1, then x* as root of (2) has a lower multiplicity. 

Wc find defiation used repeatedly first in [36] , and later modified in [33] and 
apphed in [34, 35]. 

In theory, det(A(x)) = (or maximal minors) could be used to form new 
equations. But this is neither good symbolically because the determinant is 
usually of high degree and leads to expression swell, nor numerically, as the eval- 
uation of the determinant is numerically unstable: || det(A(x)) — det(A(x))|| 



Instead of using the determinant, on a system F of N equations in n vari- 
ables, we proceed along the following three steps to form new equations: 

1. Let r = rank(A(xo)) for xq « x*. For numerical stability, we compute the 
rank via a Singular Value Decomposition (SVD) of the matrix A = A(xo). 

2. Let h e be a random vector and B G C"^*^''"'"^'' be a random matrix. 
For numerical stability, we generate random numbers on the unit circle. 

3. Let C(x) = A{x)B, notice that C = [ci,C2, . . . ,Cr+i] is an A?' x (r + 1) 
matrix with polynomial entries. 

With probability one (exceptional pairs of vectors h and matrices B belong to 
a proper algebraic subset of x C"^(''+^)) we have 



The random h and B guarantee the existence and uniqueness of the solution 
A to G(x, A) when x = x*. 

In one defiation step, we add the equations of G(x, A) instead of det{A{x)) = 
to the system F{x) = 0, adding r + 1 extra variables Ai, A2, . . ., A^+i. 

For r = n — 1, adding G(x, A) to F{x) = is equivalent to the Liapunov- 
Schmidt reduction [1] used to compute bifurcation points, see e.g. [10]. 




(2) 



||A(x)-A(x)||. 



rank(^(x*)) = r ^ corank(C'(x*)) = 1 
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Input: -F(x) = polynomial system; 
xq initial approximation for x* 



A+ := SVD(A(x,)); 
Xfc+i := Xfc - A+i^(xfe); 



r := rank(^+); 



F := 



i^(x) = 
G(x, A) = 



A := LeastSquares(G(xfe+i, A)); 
k:=k+l; Xfe:=(xfc,A); 




Output: F;Xfe_)_i. 



Figure 1: Flowchart for a modified deflation method. In practice, the test for 
quadratic convergence is implemented by a test whether the numerical rank of 
A(xfe) equals the number of columns of A{-x.k). 



3 A Bound on the Number of Deflations 

The termination of our algorithm in Figure 1 depends on the following theorem. 

Theorem 3.1 The number of deflations needed to restore the quadratic conver- 
gence of Newton's m,ethod converging to an isolated solution is strictly less than 
the multiplicity of the isolated solution. 

The answer to the question "How much less?" can be imdcrstood by looking 
at a standard basis for the ideal generated by the given polynomials in the 
system. We use standard bases to prove the termination of our algorithm, as 
explained in the next two subsections. 



4 



3.1 Standaird Bases for Local Orderings 



Let R = k[xi, ...,a;„] be the ring of polynomials in n variables with eoeffieients 
in the field k. We use the following multi-degree notation: x" = Xi^X2^ ■ ■ ■ x'^^ , 
where a = (ai, a„) is a vector of nonnegative integers. 

A multiplicative ordering < on the monoid {x" | a S Z"q} is a local ordering 
if 1 >x" for allaT^ (0,0, ...,0). 

To any weight vector w € M"o associate the weight ordering <^ by 

setting 

x" <^ x'' ^ (a,w) < {13, uj) or [ (a,w} = (/3,c^) and x" <i,^ x'^ ], (3) 

where (•, •) is the usual inner product and the lexicographic ordering, < lex, is 
used to break ties. 

In presence of a monomial ordering a polynomial 

/(x) = c^x" e R, where supp(/) = { a | ^ } is finite, (4) 

has the following attributes associated with it: 

le(/) = the leading exponent = rnaxsupp(/) 

lm(/) = the leading monomial = x'^*--'^^ 
lc(/) = the leading coefficient = Cie(/) 
lt(/) = the leading term = lc(/)lm(/) 

Let / C i? be an ideal. We call a set of polynomials S C I a. standard basis 
of / if for any / G / there is g & S such that lm((7)|lm(/). Alternatively, 5 is a 
standard basis iff the initial ideal in(7) = ({lm(/)|/ e /}) is generated by the 
leading monomials Im(S') = {lm{g)\g G S}. 

The monomials that do not belong to the initial ideal in(/) are called stan- 
dard monomials. The minimal generators of in(J) shall be called the corners of 
the staircase, the corners of the form xf for some i and a are called the endpoints 
of the staircase. 

A standard basis S is called reduced if the leading monomials of its elements 
form a minimal generating set for the initial ideal in(7) and the tail g — lt{g) 

contains only standard monomials. 

Graphically, any monomial ideal can be represented by a staircase in the 
nonnegative integer lattice Z"q. For example, let / be the ideal of i? = k[xi,X2] 
generated by 

/i = xl+xixl; 

h = xxx\^-x\\ (5) 
h = xlx2 + xixl. 
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The initial ideal depends on the ordering chosen: the staircase at the left in 
Figure 2 represents in^(/), where ui = (—1,-2). The staircase at the right in 
Figure 2 represents in^(/), where u) = (—2,-1). 



() 

C) 




() O 



x1x2 + X1X2 



( BOO 



xf + X1X2 




Figure 2: Two different staircases of the standard basis of / with respect to dif- 
ferent local orderings <(_i,_2) (at the left) and <(_2,-i) (at the right). Mono- 
mials generating in(7) are represented by black disks, while the standard mono- 
mials are shown as empty circles. 



Observe in Figure 2 that the number of standard monomials is the same for 
both orderings. This is always so, for any local ordering, as the standard mono- 
mials form a basis of the fc-linear space R(xi....,Xr^) / R{xi,....x„}-^ j where R{xi,....x„) 
is the localization of the polynomial ring R at the origin and R(xi,...,xn)^ t;he 
extension of the ideal I in this localized ring (see [5] for details). This linear 
space is of finite dimension iff the origin is an isolated solution; its dimension, 
which is the multiplicity of the origin, then equals the number of the standard 
monomials for any local ordering. 

Thanks to Mora [32], there is an algorithm for computing standard bases, 
generalized by Greuel and Pfister [16], and implemented in the computer algebra 
system Singular [17]. We will not use standard bases explicitly - except for 
theoretical purposes - but note an analytic interpretation of the local ordering 
<^: as we approach the origin along a smooth curve 



c : C ^ C" such that c{t) 



hit-^^{l + 0{t)) 



bnt-^-{l + 0{t)) 



(6) 



with uj G Z"g and (61, 6„) e C"\{0}, for every / € / the leading term lt,„(/) 
becomes dominant, i.e.: 

/(c(i)) 



\tM){c{t)) + o{f 



(7) 
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3.2 Understanding the Deflation Method 



First of all let us formulate the goal of what we would call the symbolic deflation 
process: Given a system of polynomial equations fi = 0, i = 1,2, . . .,N with 
the point x* € C" as an isolated solution of multiplicity m > 1, find a system 
Qi = 0, i = 1,2,...,N' such that x* is still an isolated solution of multiplicity 

less than m. 

The best deflation one can cook up is the one that corresponds to the max- 
imal ideal annihilating x*, i.e.: gi = Xi — x*, i = 1,2,. . . ,n. However, from a 
practical angle of numerical methods what wc actually need is an algorithm that 
would relate the deflated system to the original one in a numerically stable way 
and taking into account the fact that the isolated solution x* may be known 
only approximately. 

3.2.1 A Symbolic Defiation Method 

Here we assume that everything is exact and, therefore, w.l.o.g. we may assume 
that the isolated solution x* is the origin. 

Consider the ideal / generated by the polynomials fi of the original system. 
Wc call an ideal /' a defl,ation of / if /' D /, /' / B., and the multiplicity of the 
origin for /' is lower than that for the original ideal /. 

If the multiplicity m > 1, it means that the initial ideal in(7) does not 
contain Xi for some i. 

Proposition 3.2 Suppose m > 1 and let g be an element of a reduced standard 
basis of I with respect to a local monomial ordering <, such that lm(g) = xf, 
for some i G {!,..., n} and d > 1. Then the ideal I' = I+{dg/dxi) is a deflation 
of I. 

Proof. The derivative dg/dxi can not contain monomials > Therefore, 
/' contains / properly, since \m(dg / dxi) = xf^^ is a standard monomial for /. 
The appended generator dg/dxi still vanishes at the origin, hence, I' ^ R. □ 

Life is easy if one can compute a standard basis. However, this procedure 

is expensive symbolically and unstable numerically. Can wc find Xi and g in 
the proposition in a less straightforward way? The next lemma gives a positive 
answer. 

A linear coordinate change T : C" C" induces an automorphism of the 
polynomial ring R — C[a;i, Xn], which we call T as well: T{f )(x) = f{T{x)). 
The ideal T{I) = {T(/) | / g /} = (T(/i), . . . , T(/Ar)) represents the system 
after the change of coordinates. 

Let ^(x) be the Jacobian matrix of the system F(x) = 0, i.e.: an iV-by-n 
matrix with polynomial entries Ay (x) = dfi/dxj. The origin is singular iff 
c = corank(A(0)) > 0. Since the Jacobian matrix is rank deficient, the kernel 
of A(0) is nonzero. 
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Lemma 3.3 Take a nonzero vector X € ker A(0) C C" and let T : C" ^ C" be 
a linear coordinate transformation such that: 

n 

Ti{x) = XiXi+'^IJ,ijXj, for i = 1,2,..., n, (8) 

where [X, nonsingular matrix. 

Then di{T{I)) = {^J \ f G T{I)} is a deflation ofT{I). 

Proof ¥ov &\\i = l,2,...,N, 



J— I J J— I \ J / 

The last expression is equal to when x = 0, since A G ker A(0). 

Take any g = 6iT(/i) H h 6jvr(/jv) e T(7), where 6, G for all i. Then 

dg 9(T(/i)) g(r(/iv)) 

a — = "1 a 1 '""■'V ^ 

OXi OX\ OXi 

In view of (9), the last expression evaluates to at x = 0. Therefore, di{T{I)) 

is a proper ideal annihilating the origin. On the other hand, there is an element 
g of a reduced standard basis of T(I) with respect to a local ordering such 
that in{g) = xf with d > 1. According to the Proposition 3.2 the ideal I' = 
T{I) + {dg/dxi) is a deflation of T{I). So is di{T{I)), for it contains I'. □ 

Lemma 3.3 leads to Algorithm 1. 
Algorithm 1 G = Symbolic_Defiation(F) 

Require: F, a finite set of polynomials in R, such that the ideal {F) has mul- 
tiplicity m > 1 at the origin. 

Ensure: G, a finite set of polynomials in R, such that the ideal (G) + (F) is a 
deflation of {F). 

Compute the Jacobian ^ of F at the origin; 
Pick a nonzero vector A G ker A(0); 



3.2.2 A Numerical Deflation Method 

This method (see Algorithm 2) does boil down to the Algorithm 1 and here 
is why. 



Algorithm 2 G = Numeric_Deflation(F, xo) 



Require: F = {/i, . . . , /jv}, a finite set of polynomials in R, such that the ideal 

(F) has multiplicity m > 1 at the point x* w xq. 
Ensure: G, a finite set of polynomials in R' = R[Xi, . . . , Xr+i], where r = 

rankj4(xo), such that 

• the ideal (G) C R' has an isolated solution at the point P = (x*,A*) G 

£n,r+l. 

• the vector A* is determined uniquely; 

• the multiplicity of P is less than m. 

Compute the Jacobian matrix ^(x) of F; 

r := rankj4(xo); 

(Rl) Generate a random matrix B € ^nx{r+i).^ 
C(x) := A{x)B {N X {r + 1) matrix with polynomial entries); 
Let A = (Ai, . . . , A^+i)"^ be a vector of indeterminates; 
Consider N new polynomials ^^(x, A) = (C(x)A)j £ i?'; 

(R2) h{\) := hiXi H h hr+iXr+i — 1, where the hi are random numbers 

in C; 

G:=FU{gi,...,gN}U{h}. 



Consider a point P = (x, A) G C"'''+^ and let x = Xq. When this specializa- 
tion is performed, the values for A are determined by the following system of 
A'' + 1 linear equations: 

r ^(xo)A = 

I (h,A> = 1 ^'"^ 

where h = {hi, ... , hr+i) is a vector of random complex numbers. 

Observe how C is created: the randomization step Rl insures that the 
r + 1 columns of C(x*) are random combinations of the columns of ^(x*) and, 
therefore, corankC (x* ) — 1 with probability one. Then A is bound to live in the 
one-dimensional ker C (x* ) . The randomization step R2 makes sure one nonzero 
vector is picked out from the kernel. This proves the uniqueness of A*. 

Since C(x*)A = A{x.*)BX = 0, the vector B\ corresponds to A in Al- 
gorithm 1, provided x* is the origin. Therefore, the multiplicity drops by 
Lemma 3.3. 



4 A Symbolic-Numeric Implementation 

The method was tested and developed in Maple 9. It is implemented in PHCpack [39], 

publicly available in release 2.3. While floating-point arithmetic is used, the re- 
sult is symbolic, in the form of a new polynomial system with a well conditioned 
root. 

In our implementation we computed the full SVD. although we only need 
to know the rank of the Jacobian matrix at the current approximation for the 
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root. To speed up the algorithm, we refer to the recent techniques presented 
in [31]. 

Another issue is the exploitation of the structure in the deflated systems. 
This issue is related to automatic differentiation [19, 37]. See also [24]. 

The deflation of the system F(x) = with Jacobian matrix A{x.) leads to 
a system F^^^ (x) = whose Jacobian matrix A^^^ (x, A) has an obvious block 
structure. At the right of formula (11), the matrix A(^)(x, A) has two columns: 
the first with derivatives with respect to x and the second with derivatives with 
respect to A. 



F«(x,A) = 



F(x) = 
A{x)BX = 
hA = 1 



^W(x,A) = 



A{x) 
{^AiK))B\ A{^)B 
h 



(11) 

where the derivative of a matrix of polynomials A(x) with respect to a vector 
of unknowns x = (xi, X2, ■ ■ ■ , Xn) gives a sequence of matrices: 



dA{x.) 9A(x) 



a^(x) 



dxi 



dx2 



dx„ 



(12) 



To evaluate A^^^ (x, A) efliciently, we must first evaluate A{x) before we multiply 
with B. Otherwise, multiplying first the polynomials in A{x) with the matrix B 
will give rise to multiple occurrences of the same polynomials in A{x). Expand- 
ing A(x.)B as polynomials in x makes matters even worse. Likewise, we first 
evaluate -^A{x) before multiplying with the vector BX. 

With the exploitation of structure, the system and its Jacobian matrix can be 
evaluated about twice as fast, as our experiences on the cyclic 9-roots problem 
suggest. With 7 multipliers added to the original 9 variables, the time for 
1000 evaluations of F(^^(x, A) at a random point drops from 1.75 to 0.94 cpu 
seconds^, and the time to evaluate A(^^(x, A) 1000 times drops from 11.0 to 5.3 
cpu seconds. With 8 multipliers, the drop is even more significant: from 2.07 
to 0.92 cpu seconds for 1000 evaluations of _F'^^)(x, A), and from 12.5 to 5.1 cpu 
seconds for 1000 evaluations of A^^^ (x. A). When the rank deficiency is modest, 
the number of multipliers is close to the original number of variables and we 
can evaluate about twice as fast. 

Setting A^^' = A, B'-^^ = B, and h^-'^^ = h, it is straightforward to formally 
write down the kih. deflated system: 



fW(x,A«,...,A«) = < 



F('=-i)(x,a'^\...,A('="^)) = 

A('=-i)(x,'a«,'...,'a('=-i))BWaW = 

h(fe)A('=) = 1. 



(13) 



The benefit of exploiting the structure becomes more significant, as illus- 
trated by the evaluation of F^^\ continuing the experiment on cyclic 9-roots. 

^Execution done on a 2.4GHz Linux machine with full optimization of the code. 
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After adding 8 multipliers (corank 2), in the first deflation, we have 17 variables. 
As the corank cannot decrease, a second deflation will need 16 multipliers. The 
time to evaluate F(2)(x, A(^\ A^^^) for 1000 random values of x, X^'^\ and A^^^ 
dropped from 45.1 to 5.66 cpu seconds. 

The columns of the Jacobian matrix of F^''^ are organized in /c + 1 blocks. 
For example, for k = 2, ^(^^(x, A'^\ A^^') has in its three block columns the 
derivatives with respect to x, \^^\ and A^^^ respectively: 



(|d)5(i)AW 








h(i) 

f Mill) R(2)\(2) 










A(i)s(2)a(2) 

h(2) 



(14) 



where A = A(x), A^i) = ^(i) (x, A^^^), and means the same as in (12). No- 

tice that the first three block rows and first two block columns of A^"^^ (x, A*-^-* , A^^-* ) 
equal A(i)(x, A(^)). 



5 Applications and Numerical Results 

The implementation has been tested on seventeen examples, available 
at littp://www. math.. uic.edu/~jan/demo. html. The initial approximations 
for Newton's method were taken from the end points of solution paths defined 
by a polynomial homotopy to find all isolated solutions (see [26, 27] for a recent 
surveys). The numerical results reported in the table below are obtained with 
standard machine arithmetic. 

Observe the improved numerical conditioning in Table 1. This observation 
justifies the naming-^ of our method as a "re-conditioning" method. 

One of the interesting examples is taken from [36] and listed as "ojikaS" in 
Table 1. This system has two isolated roots: one with multiplicity two, and the 
other one has multiplicity four. Both roots need only one deflation, but at the 
double root, the rank of the Jacobian matrix is two, while the rank is one at the 
other 4-tuple root. The program produces two different deflated systems: one 
with three multipliers (for the double root) and the other with two multipliers 
(for the 4-tuple root). 

The so-called cyclic 9-roots problem is the last and largest example in Ta- 
ble 1 . This system is a widely used benchmark in the field of polynomial system 
solving, e.g.: [2, 3], [11], [12], [20], [28], with theoretical results in [21]. There 
are 5,594 (333 orbits of size 18) isolated regular cyclic 9-roots, in addition to 

the 162 isolated solutions of multiplicity four. One defiation suffices to restore 
quadratic convergence on all 162 quadruple roots of this large application. 

^We are grateful to Erich Kaltofen for this naming at the ISSAC'04 poster session. 
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System 



D 



corank(A(x)) 



Inverse Condition^ 



egl 
eg2 
eg3 
eg4 
eg5 



2-^0 
1 ^ 
1 ^ 
1 1 -> 
1 ^ 1 ^ 1 ^ 



8.3e-09 
1.2e-08 
5.6e-09 
3.0e-10 
6.4e-ll 



5.0e-01 
1.8e-01 
1.2e-01 
6.4e-02 
5.7e-03 



bakerl[22] 



1 







1.7e-08 ^ 3.8e-01 



cbmsl[38] 
cbms2[38] 
mthl91 
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3^0 
3^0 

2 ^ 



4.2e-05 
1.2e-08 

1.3e-08 



5.0e-01 
5.0e-01 

3.5e-02 



deckerl[8] 
decker2[6] 
decker3[7] 



1^1-^0 
1 1 ^ 1 ^ 

1 ^ 



3.4e-10 
4.5e-13 

4.6e-08 



2.6e-02 
6.9e-03 

2.5e-02 



ojikal[33] 
ojika2[33] 
ojika3[36] 

ojika4[35] 



1 ^ 1 ^ 
1 ^ 
1 ^ 
2^0 

1 ^ 1 ^ 



9.3e-12 
3.3e-08 
1.7e-08 
6.5e-08 
1.9e-13 



4.3e-02 
7.4e-02 
9.2C-03 
8.0e-02 
2.4e-04 



cyclic9 



1 







5.6e-10 ^ 1.8e-03 



Table 1: Numerical Results, D is the number of deflations needed to restore 

quadratic convergence. The fifth cohimn shows the decrease in the corank of 
the Jacobian matrix for all stages in the deflation. In the second to last column 
we see the change in the estimate for the inverse condition number of A(x) at 
the start of the deflation to the end of the deflation for x » x* . The last column 
lists the increase in the number of correct digits from the initial guess to the 
final approximation. 



6 Conclusion 

Our modifled deflation method works in general, is numerically stable, relatively 
simple to implement; and perhaps most importantly, a preliminary implemen- 
tation on a wide class of examples performs quite well. 

However, the doubling of the number of equations by deflation quickly leads 
to huge systems, we will search for ways to limit the number of deflations. 
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